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1 Introduction 

The theoretical understanding of the yielding transition in ather- 
mally driven disordered systems is a highly challenging problem 
and no consensus has been established even on the basic ingredi¬ 
ents that should underlie coarse grained descriptions of the non¬ 
linear rheological response of yield stress materials 

The only commonly accepted point of view is that disordered 
materials, such as glasses or soft matter, exhibit a strongly hetero¬ 
geneous dynamics when driven by an external shear. Fast particle 
rearrangements, the so-called shear transformations, take place 
in small regions while the rest of the material deforms elasti¬ 
cally^. These plastic events induce long-range elastic deforma¬ 
tions in the system leading to complex correlations of the yielding 
regions in form of plastic avalanches 

One of the major concerns remains to know whether despite 
this dynamical complexity it is sensible to describe the yielding 
dynamics within mean-field descriptions We tackle this 

key question using particle based simulations, concentrating on 
flow responses in a regime that is relevant for many rheological 
setups We find that at high enough shear rates and/or 

small enough system sizes we recover a dynamics well described 
by mean-field considerations^^, similar to near mean-field critical 
points in equilibrium phase transitions 

An important feature that found its way into several theo¬ 
ries, e.g. the Soft Glassy Rheology (SGR) and the Shear Trans¬ 
formation Zone Theory (STZ), is to describe yielding in a gen¬ 
eralized thermodynamic description with an effective tempera¬ 
ture . But the analogy with thermally activated events proposed 
by the SGR theory has recently been questioned for athermal rhe- 
ology^’^. Also it has been shown that the thermodynamic inter¬ 
pretation of the effective temperature is problematic in the ather¬ 
mal regime and that the STZ theory is not able to predict the 
non-linear response expected in the small driving limit 
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Fig. 1 Schematic view of the coarse grained picture: The particle 
simulation box is devided into smaller mesoscopic parts that can hold 
exactly one plastic event. The plastic activity in the surrounding of a 
mesoscopic region (see graphical interpretation on the right) leads to 
stress fluctuations around the mean value of the stress (dashed line in 
the graph on the right), which is controlled by the external forcing. 


In this work we show that the dynamically created noise in 
sheared amorphous systems can be encompassed through a nor¬ 
mal diffusion in the local stresses (for a schematic view see 
Fig. 1). In agreement with former works we consider that 
local stress diffusion is acting as noise in the tilt of the local po¬ 
tential energy landscape with an amplitude proportional to the 
plastic activity, an observable that can be measured experimen- 

tally2h24_ 

We analyse within molecular djmamics simulations the cou¬ 
pling strength between the mechanical noise and the plastic ac¬ 
tivity. Within our simulations we successfully relate this cou¬ 
pling strength, a dimensionless and density independent quan¬ 
tity, in a consistent manner to the flow response without any 
further parameters. This finding leads us to reconsider a spec¬ 
ulative mean-field scenario put forward a long time ago in the 
so-called Hebraud-Lequeux (HL) model that was at the basis 
of important further developments 1^,22,34^ providing so far one 
of the best self-consistent mean-held description of mechanical 
noise in athermally sheared disordered systems. 

2 Mean-field approach 

Within the framework of the HL-model we can establish a link be¬ 
tween the dynamical yield stress and the prefactor of a Herschel 
Bulkley type power law ht of the rheological curves of athermally 
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sheared yield stress materials. This relation should, according 
to this theory, solely be determined by the coupling strength be¬ 
tween mechanical noise and the rate of activity. In the following 
we review briefly the main assumptions and results of the mean- 
field approach^. The model describes the state of a soft glassy ma¬ 
terial via the probability density of local shear stresses a = oVj, 
in regions of mesoscopic size IT while the material is sheared at 
rate y. The time evolution of is given by 

d,g^{o,t)= -GQy{t)da&’-\ei\o\-ac)S^ 

+r(05(cT)+DHLW<3^^ (1) 

where 0 (x) and 5 (x) denote respectively the usual Heaviside and 
delta-distributions. The first term on the right hand side pro¬ 
portional to the stress gradient of the probability density da^ 
accounts for the linear elastic response. The following term de¬ 
scribes the loss in the probability density due to local yielding of 
overstressed regions above a critical stress Oc at a rate given by 
l/x. It has been argued that the phenomenological parameter f 
can be interpreted as the duration of a plastic event in the low 
shear rate limit The corresponding gain term is given in the 
third expression on the right hand side, where the stress is set to 
zero after a yielding event with a rate given by the plastic activity 
rate 

r{t) = - [ ‘^{o,t)da. ( 2 ) 

The last term encompasses the stress changes created through 
other yielding events in a mean-field manner, assuming that this 
mechanical noise can be approximated through a normal diffu¬ 
sion of the mesoscopic stresses. To describe this noise in an self- 
consistent manner, the HL approach proposes that its amplitude 
should be related to the rate of plastic activity through a dimen¬ 
sionless coupling constant a 

DHL(f) = aoMt). (3) 

This last relation introduces a non-linearity into Eq. (1), since 
the rate of plastic activity itself depends on the density probability 
of the mesoscopic stresses. It is this coupling that renders the 
problem non-trivial and yields interesting results regarding the 
behaviour of macroscopic quantities. 

This model is known to exhibit a unique stationary state for 
a finite shear rate in the large time limit, where the probability 
density for the mesoscopic stresses becomes time independent. 
To determine the time averaged macroscopic stress in the steady 
state one averages over the mesoscopic stresses weighted by the 
corresponding steady state probability density 

{a) = j (T^{a)do. ( 4 ) 

Using appropriate units we can write the equations in dimension¬ 
less quantities, expressing stress related values in units of the lo¬ 
cal yield stress etc, time quantities in units of t, the shear rate in 
units of ctc/(Got) and the stress diffusion coefficient in units of 
CT^/f, leaving only two independent dimensionless model param¬ 
eters that determine the flow behaviour, namely the dimension- 



Fig. 2 Flow curves: Rescaled macroscopic shear stress o-j^/etc as a 
function of the rescaled applied shear rate y(Got:o/oc) for different 
densities (see text). The dashed line is a fit with the Herschel-Bulkley 
form (Tej, = (jy -t Ay'/z. jhe inset shows the same data before rescaling, 
dashed lines display a guide to the eye with the same Herschel-Bulkley 
fitting form. 

less shear rate and coupling constant a. 

The rheological results in the small shear rate limit for this 
model are well studied For small enough coupling 

strength a < 1 /2 the HL model predicts a Herschel-Bulkley flow 
behaviour of exponent 1/2, (a) pa cr^+Ay*/^, for the time aver¬ 
aged macroscopic stress in the steady state (rr) = / aS^(c>)da, 
with the two constants Oy (the dynamical yield stress) and A (the 
prefactor). These macroscopic constants of the model can only 
depend on the last free control parameter, namely the coupling 
constant a. This means that we will obtain a universal curve* 
for the rescaled quantities A/^GqTGc versus Gy/G c parametrized 
through a (see Fig. 3). 

The physical interpretation of this coupling constant could be a 
mechanical fragility^ that would depend on the details of the mi¬ 
croscopic interactions and thus should be material dependent. In 
a later work by Boequet et ah, the authors use a spatial approach 
to the problem that allows for the derivation of an expression 
for the coupling constant a if one assumes a decorrelation of the 
plastic event dynamics. Within this approximation it is possible 
to express a as a function of the specific form of the elastic stress 
propagator This suggests that there should be classes of ma¬ 
terials with similar elastic responses that not only share the same 
non-linear exponent in the flow curve, but also comparable rela¬ 
tions between the dynamical yield stress and the prefactor in the 
flow curve. 

In the following we aim at probing not only the above con¬ 
stant relation using microscopic dynamic simulations, but also to 
test the the underlying assumptions, most importantly the self- 
consistent description of the mechanical noise induced through 
the plastic activity. To test the above theory we need to measure 
all involved parameters that appear in the HL description. 

Quantities and parameters like the local stress and the shear 


* For the analytical derivation of this constant relation see the Appendix A 
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Fig. 3 HL predictions vs. MD simulations: (a) The points are the rescaled 
Herschel Bulkley prefactor A as a function of the rescaled dynamical 
yield stress ay for four different densities, error bars are estimated from 
the measurements of the different quantities. The full line corresponds 
to the HL prediction, (b) This inset shows the data from the simulation 
fits before the rescaling, (c) This inset shows a zoom into the main 
panel, to display better the proximity to the theoretical curve obtained 
from a very simplified picture. 

modulus Go are well defined quantities in the microscopic sim¬ 
ulations and rather easy to measure. Also, if stress fluctuations 
turn out to be diffusive, £>hl has a well defined meaning. How¬ 
ever, other HL parameters like the local yield stress etc, the phe¬ 
nomenological parameter f and the rate of plastic activity r(r) are 
rather difficult to interpret within the microscopic picture. In the 
following we describe our attempt to relate the different parame¬ 
ters to measurable quantities in the microscopic dynamics. 

3 Microscopic modei 

We have investigated a generic two-dimensional (2D) model of 
a glass, consisting of a mixture of A and B particles, with = 
10400 and Nb = 5600, interacting via a Lennard-Jones potential 
^apir) = 4e«/j[(cT„^/r)'2 - {Oap/rf] with a.jS = A,B and r be¬ 
ing the distance between two particles. The parameters £^(, 4 , 0^4 
and define the units of energy, length and mass; the unit of 
time is given by Tq = oAa V'('”A/eAA)- We set Eaa = 1-0, Eab = 1-5, 
£bb = 0.5, Gaa = 1-0, oab = 0.8 and cjbb = 0.88 and = mb = I- 
This choice is known to prevent crystallization in 2D at low tem¬ 
perature'^^. The potential is truncated at r = rc = 2.5 for compu¬ 
tational convenience and periodic boundary conditions are used. 
The equations of motion are integrated using the velocity Ver- 
let algorithm with a time step St = 0.005. The athermal limit is 
achieved by thermostating the system at zero temperature via a 
Langevin thermostat with a damping coefficient ^ = 1 which 
corresponds to a strongly overdamped condition for the dynam¬ 
ics'^. This model, which has been widely studied in different 
versions is usually considered appropriate for colloidal 

and other soft glasses. This choice is motived by the purpose of 
investigating the general aspects of the rheology of athermal yield 
stress materials. 

To investigate different athermal flow responses, we explore 
states with different number density p = (A (4 +Nb)/V by chang¬ 


ing the volume V of the system. Glassy configurations were pre¬ 
pared by quenching equilibrated configurations at T = 1 to zero 
temperature with a fast cooling rate. Simple shear is set at a rate 
y by deforming the box dimensions and remapping the particle 
positions. The quenching protocol has virtually no effect since 
we focus on the steady state shear (total imposed deformation 
Ay >20%). 

In order to characterize the plastic activity of the system we 
consider the quantity^. For a given particle i, is defined 
as the minimum over all possible linear deformation tensors e of: 

St) = ^ [r/y(? + 8t) - (I + e) • r/y(?)] ^ (5) 

J 

where the index j runs over all the neighbors of the reference 
particle i and I is the identity matrix. We set the time lag to 
St = 4. This value is a compromise between having a good sig¬ 
nal, i.e. large irreversible displacements, and being able to resolve 
individual plastic events. 

3.1 Macroscopic flow curve 

In Fig. 2 we show the dependence of the macroscopic shear stress 
(Jxy on the applied shear rate y. The flow curves are well described 
by the Herschel-Bulkley (HB) law, Oxy = Oy +Af', with an expo¬ 
nent n Ks 0.55 (not shown here) which seems not to depend on the 
density. Fixing the exponent n to the value 0.5, the one predicted 
by the HL model in the case of a < 1 /2, gives indistinguishable 
results. Although other works on sheared disordered material in 
two and three dimensions report similar values for the flow curve 
exponent recent works in the literature seem to suggest that 
a proper finite size scaling analysis close to the yielding transi¬ 
tion reveals different critical dynamical exponentsWe 
would like to insist here on the fact, that our study does not aim 
at measuring the critical exponents of the transition; instead we 
rather test the consistency of the assumptions made in the HL ap¬ 
proach in a parameter regime that fits well the model predictions. 

3.2 Size of an elementary plastic region 

Here we describe the procedure we followed to convert the micro¬ 
scopic simulations into a mesoscopic description. First, we denote 
a particle i as active, i.e. performing a plastic rearrangement, at a 
given time t if the corresponding D^j^^{i,t) is larger than a thresh¬ 
old value that we fix equal to 0.1. In Fig. 4(a) we compare a typi¬ 
cal stress-strain curve and the corresponding evolution of the total 
number Npi of active particles. The correlation between the stress 
drops and the peaks in Npj suggests that the actual definition is 
reasonable. Then a high-resolution discretization of the system is 
performed by dividing it into kxk square blocks of length wq = 2. 
A small block is considered as active if it contains at least one 
active particle. 

The mean size, i.e. the mean linear extension, of a plastic 
event (Ipi) can be estimated by a cluster analysis of the spatial 
arrangement of the active blocks in the configurations explored 
by the system. We employed a modified version of the Hoshen- 
Kopelman algorithm’ in order to account for periodic boundary 
conditions. If we assume (Ipi) = (Api)^/^ where Api is the area of 
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a plastic cluster, we obtain {Ipi) ss 6 , with no relevant dependence 
on the density and shear rate (in the range 7 ^ 10 The value 
is in accord with previous works reporting plastic regions with 
a size of a few particle diameters Furthermore, this agree¬ 
ment justifies the criterion we employed to define active particles. 
Indeed, if we improve the resolution in the analysis of the plas¬ 
tic activity by decreasing the threshold on by a factor 10 , 
while the number of active particles increase by a factor 4, the 
mean extension of a plastic event {Ipi) is reduced by half, suggest¬ 
ing that single particle rearrangements are erroneously taken into 
account. 

Next we implement the coarse-graining of microscopic simula¬ 
tions on the scale of individual plastic events. The simulation box 
is divided into M x M square blocks with M chosen in order to 
have IV ^ {Ipi)- The local shear stress o’™ of a block m is defined 
as: 


_m __ ^ 




dV{r‘^) r'Jry 
dr^J r^J 


i£m j=l 


( 6 ) 


where vi^x and vij are the v and y components of the velocity of 
the particle i, r‘j is the distance between the particles i and j and 
the summation of i is performed over the particles in the block. 
The macroscopic stress tensor o^ is obtained by the summation 
of o™ over all the blocks. 



Fig. 4 Activity and stress diffusion: (a) The full line represent a part of 
the macroscopic stress-strain curve in the steady state regime for a 
density p = 1.149 and a shear rate 7= lO^®. The average slope of the 
elastic parts on the curve yields Go- Open symbols: corresponding 
evolution of the number of active particles Npi in the system, (b) Coarse 
mean-square stress difference ((Ao^j^jf)) as a function of time for 
different values of the density. Note the linear behaviour at short times. 
Inset: stress diffusion coefficient Da = {{Aaf^)'^{t))/t as a function of 
density. 


3.3 Stress diffusion and duration of a plastic event 

The above introduced mean-field model assumes local stress fluc¬ 
tuations, obeying a normal diffusion process in stress space. To 


test this idea we define the coarse mean-square stress difference 
as: 

((A(7^;)2(f)) = ((cT™(fo + f)-a:;(fo)-Gon)') (7) 

where C7™(f) is the stress in a given block at time t and the last 
term in Eq. (7) accounts for the stress increase due to the elas¬ 
tic deformation of the system, being Gq the macroscopic shear 
modulus. In Fig. 4(b), the short time beha’vior of ((Ao’"^)^(f)) 
is shown for a finite shear rate 7 = 10 ^®, approaching the qua¬ 
sistatic limit. We observe that at short times, ((A( 7 (p^(f)) in¬ 
creases linearly with time. We define the stress diffusion coef¬ 
ficient as Dr, = {{A<j"y}'^{t))/t. In the inset of Fig. 4(b) we show 
the dependence of Da on the density of the system. 

To define a duration of a plastic event we analyse the two- 
time autocorrelation function of the Df^-^ quantity for active 
particles. In Fig. 5(a) we show this correlation function Cp = 
(^L(fo)£»L«(fo + 0)/((£’L’«(fo))^) as a function of time. We ob¬ 
serve that Cp decays exponentially with a characteristic time Tp, 
that depends weakly on density (see the inset of Fig. 5(a)) and 
that is close to the damping time ^ * of the Langevin ther¬ 

mostat. We choose to interpret this decorrelation time as the typ¬ 
ical duration of a plastic event entering the HL model description. 

3.4 Local yield stress 

In this section, we present a method which allows us to calculate 
a local critical stress, i.e., the stress limit before a plastic rear¬ 
rangement occurs locally. For this purpose, we adopt the “frozen 
matrix" method The system is frozen except for a target re¬ 
gion, i.e., a mesoscopic block, and it’s subjected to a simple shear 
deformation, with a quasi-static protocol. The frozen region can 
only deform affinely whereas the target block is allowed to relax 
non-affinely. For small deformations the target region behaves 
elastically and the stress increases linearly, with a slope controlled 
by the local shear modulus^®. As the strain increases, the elastic 
beha’vior goes on until the local yield stress is reached and a plas¬ 
tic rearrangement takes places. This is indicated by a stress drop 
in the stress-strain curve. For a given block, we define as the local 
yield stress cr™ the value of the local stress at the first maximum. 

In Fig. 5(b) we show the distribution of a'" for the different 
values of the density. First, the fact the local yield stress is dis¬ 
tributed is in clear contrast with the HL assumption of a unique 
critical stress Oc- The exposed distributions are well described by 
a simple Gaussian forms. The mean values (< 7 "') and the variance 
{So"') of the distributions are shown in the inset of the figure. As 
the density increases, (o’™) increases, due to the enhancement of 
the repulsive interactions between particles. 

With the frozen matrix method we estimated also the stress 
release following a local yielding event. We observe that the plas¬ 
tic event only partially relaxes the accumulated stress in contrast 
with the assumption of the HL model of a complete relaxation. 
The fraction of the relaxed stress seems not to depend on density 
being (Ao™)/(o™) ^ 0.2 (see inset of Fig. 5). 

In a former work^° it was shown that the mean value of the 
shear modulus obtained with the frozen matrix method depends 
on the target region size W and that it converges, from higher 
values, to the macroscopic modulus as W increases. This is due to 
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Fig. 5 Duration of events and local yield stress: (a) autocorrelation 
function for active particles Cp = (£>L(^o)£>L>,(fo + 0>/((£>L>,(fo))^) as a 
function of time. Lines are fitting curves exp(-r/T,,). Inset: dependence 
of Tp on the density, (b) Distribution of the local yield stress p{a^’) for 
different values of the density. Full lines are Gaussian distributions fits. 
Left inset: Mean value (cr"’) and variance {Set"’) as a function of density. 
Dashed lines are power-law guides. Right inset: Absolute local stress 
drop (Act"') and the relative local stress drop (Act'")/(ct™) as a function 
of density. Dashed lines are power-law and constant value guides to the 
eye. Symbol code is the same as in panel (a). 


the frozen environment which reduces the non-affine motion of 
the particles in the target region. For W ^ 6 the discrepancy is up 
to 50%, a significant error. We are aware that the estimate of (ct'") 
may be affected by this effect. Also we have to be careful in the 
interpretation of the result, because by using the frozen matrix 
method we measure the yield stress distribution, obtained from 
strained configurations. It would be interesting to try to infer the 
inherent yield stress distribution from our measurements®, but 
this is beyond the scope of this study and left for further investi¬ 
gations. 

4 Robustness of the HL model 

In the following we aim at probing the robustness of the HL model 
by comparing our data from the microscopic simulations to the 
various assumptions and predictions of the mean-field descrip¬ 
tion. 

Let us first recall some of the basic assumptions made to write 
the evolution Eq. (1) for the probability distribution of the meso¬ 
scopic stress: First (a) the yielding happens at a constant homoge¬ 
neous yield stress Gc- Then (b) the local response to surrounding 
plastic events is encompassed through local diffusive stress fluc¬ 
tuations with a well defined diffusion coefficient DHiit)- And the 
third assumption (c) concerns the relaxation of the local stress to 
zero, once a site yields, leading to a typical stress jump Aghl ~ <7c- 
The HL model assumptions, aiming for a self-consistent descrip¬ 


tion of the mechanical noise, thus predict a diffusion coefficient 
proportional to the square of this typical local stress jump Aghl 
divided by a typical time scale, given by the inverse plasticity rate 

no- 

(a) Our data analysis reveals that a homogeneous yield stress is 
of course not verified in a disordered system, where one expects 
a distribution of yield stresses (see Fig. 5b). However, it has been 
recently shown that the existence of a local yield stress distribu¬ 
tion in the HL dynamics does not strongly alter the predictions for 
the flow behaviour®. 

(b) We tested as well the second assumption of a normal diffu¬ 

sion of the mesoscopic stresses, and we find that within the range 
of shear rates that we consider, the measurement of the average 
mean-square stress differences indeed allows for the determina¬ 
tion of a well-defined diffusion coefficient. The HL coupling con¬ 
stant = Da/iGpT) (see Eq. (3)) with Da from the results 

in Fig. 4(b) and F from the cluster analysis turns out to be very 
small, with of order 10^^ for all considered densities. 

(c) However, our study suggests that the partial relaxation of 
the local stresses after yielding, {Ag"') /{g"‘) ^ 0.2 (instead of the 
assumed total relaxation) introduces an important corrective fac¬ 
tor for the coupling constant. Since the typical stress jump in the 
diffusion process is now reduced, we obtain for the stress diffu¬ 
sion coefficient an altered expression 

OCT = «micro(A(T'")2r. (8) 

We report the values obtained for this new definition of the cou¬ 
pling constant in Table 1. It turns out that these values are in good 
approximation density independent and yield an average value of 
about a pa 0.26. Thus the coupling constant displays a reasonable 
value smaller than 0.5 for which the HL model predicts® a well- 
established Herschel Bulkley regime in the flow curve with an 
exponent of 1/2 . 

And indeed our flow curves from the microscopic simulations 
can be well fitted with an Herschel Bulkley expression Gxy = 
CTy-LAy*/^ (see inset of Fig. 2). In Fig. 2 we plot our data in 
dimensionless units, rescaled as suggested by the HL model, us¬ 
ing dimensionless quantities GxyjGc and y(Got/etc). The shear 
modulus Go can be easily accessed in simulations as the slope 
of the stress-strain curve in the elastic regime. We approximate 
the critical stress Gc as the mean value (cr™) of the local yield 
stress obtained with the frozen matrix method and we associate 
the phenomenological parameter t with the duration of a plastic 
event in the low shear rate limit®. Hence, we will approximate t 
by the value of Tp obtained through the measurement of the 
two-time autocorrelation function. We observe a collapse of the 
flow curves for different densities. This collapse is not perfect, 
but regarding the large error bars introduced by our methods and 
estimations it appears still convincing and suggest a very generic 
flow behaviour. 

One of the strongest prediction of the HL model is that both 
the prefactor A as well as the dynamical yield stress Gy are de¬ 
termined solely by the specific value of the coupling constant be¬ 
tween diffusion and plastic activity. In Fig. 3 we compare the 
dimensionless yield stress ct,./(< 7™) and prefactor A/(GoCTcT)*/^, 
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density 

^m/cra/10 


1.149 

2.81±0.21 'I 

3.21±0.12 'I 

1.189 

2.47 ±0.27 ^ 

3.30±0.19 o ..o , M 


> 2.57 ±0.26 

> 3.28±0.16 

1.253 

2.27 ±0.24 

3.28±0.18 

1.322 

2.74±0.30 ^ 

3.34±0.18 ^ 


Table 1 Comparison of different coupling constant measurements as explained in the text. 


obtained from a Herschel Bulkley fit of the flow curves, with the 
analytically obtained parametric curve. We observe that all the 
data belonging to different density values collapse as expected 
by the previous collapse of the data roughly onto a single point 
close to the theoretical curve, that lies within the estimated er¬ 
ror. This result points to a universal determination of the flow 
curves, a priori strongly sensitive to the density (note the density 
dependence of A and cr,. in the inset of Fig. 3), through one single 
density independent parameter, namely the coupling constant a. 
As suggested by the KEP model we expect this parameter to be 
only dependent on the specific form of the elastic propagator^. 

If we assume the HL parametric relation between prefactor 
A (a) and dynamical yield stress ay (a) to hold, it is possible to 
estimate the coupling constant a in an alternative way through 
the macroscopic measure of the flow curve, in the following re¬ 
ferred to as af[ow Despite all the rough approximations we had to 
make, that tend to introduce large error bars on the data, we find 
the comparison between the coupling constant a,„icm with the al¬ 
ternative measurement using the expression for the rescaled yield 
stress^ ay/ac{Sc) rather convincing (see table 1). Both measure¬ 
ments suggest a density independent result with a quite small 
relative error of approximately 20%. 

Altogether our data suggests that revisiting the rule of setting 
the stress to zero after a yield event in the HL model equations, by 
changing the gain term in the evolution Eq. (1) together with the 
introduction of a more realistic yield stress distribution, seems to 
be a promising route to reach a more realistic mean-field model¬ 
ing of athermally sheared amorphous systems. 

5 Conclusion 

In this study we aimed at testing some of the most basic assump¬ 
tions and predictions of mean-field modeling for the rheology of 
athermally sheared amorphous systems. In conclusion we obtain 
a consistent picture of how to model correctly the mechanical 
noise in the regime of large enough driving rates (far from the 
true critical point We find that we can incorporate the noise 
into a normal diffusion of local stresses with a noise amplitude 
solely governed by the rate of plastic activity as proposed by the 
Hebraud-Lequeux (HL) model We not only confirm this phys¬ 
ical picture using molecular dynamics simulations, but we also 
show that the coupling strength between diffusion and the rate of 
plastic activity, a dimensionless and density independent quantity, 
seems to determine the specific form of the rheological response. 
Our data analysis suggest some important modifications in the 


t given in Appendix A 


original version of HL model equations such as the partial relax¬ 
ation of the local stress after a yielding event and the introduction 
of a yield stress distribution obtained from the microscopic simu¬ 
lations. In a future work we plan to test the coherence of such a 
modified model with our microscopic approach. Further it would 
be highly desirable to test the degree of generality of the results 
on other model systems and experimental setups, that have ac¬ 
cess to the measure of local plastic activity 
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Appendix A 

The parametric relation between the yield stress and the HB 
prefactor 

In this section, we compute the relation between the dynamic 
yield stress ay and the prefactor A in the relation (a) = ay TAy*/^ 
derived from the HL model (as given in Fig. 3). We use the 
method developed in^®: the probability density function solving 
the (stationary) HL equation (2) of the paper is expanded in the 
following way 


^{a) 

= Q°(a) + f/^Q^(a) + --- , 

for a in [-cTc,o'c] 

(9) 

^(a) 


, for ±CT in [etc, ±°°[ 

(10) 


This ansatz has been proved to be correct in the case where 
a < ct^/2 (again see®®). Moreover, term-by-term integration is 
allowed which means that 

(a) = j aB^{a)da = aQ°{a)da^ 

Oy 

+ (^l aQ^{a)da'^f^^ + --- ( 11 ) 

'-V-' 

A 
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where the dots are terms of higher order than 7 */^. Now all that 
is left to do is to identify and 2*. By plugging (9) in (1) we 
obtain the following equation on : 


-dd^Q^ + rGQdaQ^ = iSQ 
Q^\±Oc)=0 
,rk2°(ct)dcT= 1 


( 12 ) 


Note that in this equation d is an unknown coefficient used to 
enforce the integral condition on 2 *'; it is physically related to the 
diffusion coefficient by tD ~ dy. 

This system can easily be integrated and we find 2® to be 


e”(tT) = 


1 


-dcGgr 

Orll—e d 


(TGqT —OcTGq 

[ — e 3 

(cr-(Tc)GoT 

1 — e 3 


-CTc < O’ < 0 


0 < (T < etc 

and d is selected so that the following equation holds true: 

ct 2 OcGqX \ 2d J 


(13) 


(14) 


Note that this equation (in d) has a unique solution d{a) if and 
only if 0 < cc/a^ < 1/2. However, in view of (13), it is easier 
to express everything in terms of the parameter d instead of a: 
then the limit a/cr/ —> 1/2 is equivalent to d —> +00 and a —> 0 is 
equivalent to d —> 0 . 

Now using the formula for computing cTy we obtain by integra¬ 
tion. 


2.d J GcGqT 


(15) 


To obtain the prefactor A we must compute 2* which necessi¬ 
tates the computation of Using the continuity of d(,^ at ±ac, 
one can find that R]_ (z) are functions satisfying 


{-dd^R^±+Rl=0 
d,R\{0) = d„Q^\o,) = - 1 %^ 

dadl—e 3 ) 

-(TetGo 

d,Rl{0) = d^Q°{-a,) = 

daJl—e d ) 


(16) 


Again, it is easy to solve this system: 

tGo 


<(z) = 


\/dac(i~e 3 “) 


-eVd 


-gctGp 

1 xGqc 3 ^ 

(z) = ^ 

\/dac{l-e 3 ) 


(17) 


(18) 


Now we can write down the equations satisfied by 2* using the 
continuity of at <7 = ±etc: 

{-ddlQ^+GoxdaQ^ =d{l.da^dlQ^) 


Q\ac)=R\(Qi)= 

VdGc{l—e d ) 

—gcTGp 

Q\-cJc)=Rl{Q) = 

'/dGc{\—e d ) 

t/_\2'(g)dc7 = 0 


(19) 


In this system d is an unknown coefficient to be simultaneously 
computed with 2* ■ The computation of 2* is tedious but straight¬ 
forward. The expression of are quite lengthy but can be 
checked out on a symbolic computation program: if — etc < g < 0 
then 


o'(3)- 

acVd[l-e^i-j 


e 3 


dGoT , , ?alfj 

T -H5^(g + gc)e3 

d^^Gr \l — e d 


( 20 ) 


e d ~ e d 


cj^GoT / d "s/d 


d 

GqT ( 


dCr 0 } / da. 


+ - 




and if 0 < tr < (Tc 

Q\^)= - 


-f 


-b 


dG^T 

d^ac{l-e~^ j 


_GgT 

e d — e d 


(TcGo'T I d ^/d 
dOr g} 


^ 0 ^ ( 




d 

Co 
1 — e 3 




do. 


T(fT-(T, 


( 21 ) 


The parameter d is selected to enforce the vanishing integral con¬ 
dition which amounts to taking 


d = -- 


dV2 e-^ -e--^ + 

go 


( 22 ) 


Finally we can compute the prefactor A which is equal to 
/<t 2* (c!')d(7. All in all, we obtain (let us note u = {xGqOc)/ d): 


— = - coth (u) — - 
CT,. 2 u 


A 1 e“ (MCOsh(2M)+2sinh(2M)+ m) 


\/OcT:Gq 2y/u M (sinh (2m) — m) 

1 6m2 coth (2m) 

2-y/M M (sinh (2m) — u) 


(23) 


(24) 
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